/* This file uses GSL - GNU Scientific Library
You can get it from: 
http://www.gnu.org/software/gsl/
*/

#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <iostream>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h> 
#include <string.h>
     
using namespace std;

 int main( int argc, char **argv )
{
       ifstream data1;
       data1.open("data",ios::in);

       if(data1==NULL){
	        fprintf( stderr, "wrong training file! \n");
      		exit( -1 );
       }

       int num_inputs;
       data1 >> num_inputs;

       int counter = 0;
       while(!data1.eof()){
	  double temp;
          data1 >> temp;
          if(data1.eof())
            break;
//          cout << temp << " ";
	  counter ++;
//          if(counter%(num_inputs+1) == 0)
//		cout << endl;

       }

       data1.close();

//       cout << "reading from the file is finished" << endl;
//       cout << "counter: " << counter << endl;

       ifstream data;
       data.open("data",ios::in);
       double alaki;
       data >> alaki;
       
       int size_data = counter/(num_inputs+1);
 //      cout << "No. data points: " << size_data << endl;

       double x[num_inputs*(size_data)];
       double y[size_data];

       for(int j=0;j<size_data;j++){
		for(int i=0;i<num_inputs;i++){
			double temp;
               		data >> temp;
			x[j*num_inputs+i] = temp;
		}
                double temp;
	        data >> temp;
          	y[j] = temp;
	}


/*       for(int i=0;i<num_inputs*size_data;i++){
        	cout << x[i] << " ";
                if((i+1)%num_inputs ==0)
 			cout << y[i/num_inputs] << endl;
       }
*/   
       double c[num_inputs];
       for(int i=0;i<num_inputs;i++) c[i] = 0;
       double xx[num_inputs*num_inputs];
       for(int i=0;i<num_inputs*num_inputs;i++) xx[i] = 0;
       double xxx[num_inputs*size_data];
       for(int i=0;i<num_inputs*size_data;i++) xxx[i] = 0;


       gsl_matrix_view A = gsl_matrix_view_array(x, size_data, num_inputs);
       gsl_matrix_view AA = gsl_matrix_view_array(xx, num_inputs, num_inputs);
       gsl_matrix_view AAA = gsl_matrix_view_array(xxx, num_inputs, size_data);
       gsl_matrix_view B = gsl_matrix_view_array(y, size_data, 1);
       gsl_matrix_view C = gsl_matrix_view_array(c, num_inputs, 1);
     
       /* Compute AA = X'X */
     
       gsl_blas_dgemm (CblasTrans, CblasNoTrans,
                       1.0, &A.matrix, &A.matrix,
                       0.0, &AA.matrix);


      /* calculate the inverse of AA */
      gsl_matrix* inverse = gsl_matrix_alloc(num_inputs, num_inputs);
      gsl_permutation* perm = gsl_permutation_alloc(num_inputs);
      int s=0;
      gsl_linalg_LU_decomp(&AA.matrix, perm, &s);
      gsl_linalg_LU_invert(&AA.matrix, perm, inverse);


     
     /* calculate AAA = inv(x'x)x' */
     gsl_blas_dgemm (CblasNoTrans, CblasTrans,
                       1.0, inverse, &A.matrix,
                       0.0, &AAA.matrix);

     /*calcute c = inv(x'x)x'y */
     gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                       1.0, &AAA.matrix, &B.matrix,
                       0.0, &C.matrix);
     
/*     for(int i=0;i<num_inputs;i++)
       cout << c[i] << " ";
     cout << endl;

     for(int i=0;i<size_data;i++){
	double sum = 0;
	for(int j=0;j<num_inputs;j++){
        	sum += c[j]*x[i*num_inputs+j];
	}
	cout << sum << ":" << y[i] << endl;
     }
*/
      ofstream fid; 
      fid.open("weights.txt", ios::out);
      for(int i=0;i<num_inputs;i++)
	 fid << c[i] << "\t";
      fid.close();
      
      char command[1024];
      strcpy(command,"rm -f weights-");
      strcat(command,argv[1]);
      strcat(command,".txt");

      cout << command << endl;
      if(system(command) != 0){
	fprintf( stderr, "The learning part is failed!\n");
	exit( -1 );
      }

      strcpy(command,"cp weights.txt weights-");
      strcat(command,argv[1]);
      strcat(command,".txt");
      cout << command << endl;
     if(system(command) != 0){
	fprintf( stderr, "The learning part is failed!\n");
	exit( -1 );		
      }

      return 0;  
     }

